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A new computational algorithm, the discrete singular convolution (DSC), is introduced for com- 
putational electromagnetics. The basic philosophy behind the DSC algorithm for the approximation 
of functions and their derivatives is studied. Approximations to the delta distribution are constructed 
as either bandlimited reproducing kernels or approximate reproducing kernels. A systematic proce- 
dure is proposed to handle a number of boundary conditions which occur in practical applications. 
The unified features of the DSC algorithm for solving differential equations are explored from the 
point of view of the method of weighted residuals. It is demonstrated that different methods of 
implementation for the present algorithm, such as global, local, Galerkin, collocation, and finite 
difference, can be deduced from a single starting point. Both the computational bandwidth and the 
accuracy of the DSC algorithm are shown to be controllable. Three example problems are employed 
to illustrate the usefulness, test the accuracy and explore the limitation of the DSC algorithm. A 
Galerkin-induced collocation approach is used for a waveguide analysis in both regular and irregular 
domains and for electrostatic field estimation via potential functions. Electromagnetic wave prop- 
agation in three spatial dimensions is integrated by using a generalized finite difference approach, 
which becomes a global-finite difference scheme at certain limit of DSC parameters. Numerical 
experiments indicate that the proposed algorithm is a promising approach for solving problems in 
electromagnetics. 

Keywords: Discrete singular convolution; Maxwell equation; computational electromagnetics; 
waveguide; potential function; wave propagation. 



*Tel; (65) 874 6589; Fax: (65) 774 6756; E-mail: cscweigw@nus.edu.sg 



1 



I. INTRODUCTION 



Recently, computational electromagnetics (CEM) has emerged as a distinct scientific discipline for the study and 
understanding of a wide variety of electrical and electronic engineerings problems [p]-|l0| . As a natural extension to 
the analytical approach to the Maxwell equation, the CEM is based on numerically solving the governing equations 
in either the partial differential form or in the integral equation form. The complexities of physics and of the 
field geometry are no longer the limiting factors to CEM as they are to the analytical approach. With the advent 
of high-performance digital computers, CEM is emerging as a powerful approach for solving practical problems in 
electromagnetics. In fact, the dramatic progress in solving the Maxwell equation made in the last few decades l^-[ic|] 
has opened up a new research frontier in electromagnetics, plasmadynamics, optic engineering, and an interphase 
between electrodynamics and quantum dynamics. As a computational discipline, the success of the CEM is vitally 
dependent on the virtues of computational algorithms, such as numerical accuracy, stability and efhciency. These in 
turn depend on grid methods and numerical schemes for the solution of the Maxwell equation. 

A variety of computational techniques have been used for CEM, including wavelet analysis integer lattice 

gas automata [0 , hierarchical tangential vector finite elements Q] , Nedelec tetrahedral element method Q and other 
approaches |^,^. Typically, grid methods used in CEM are either global pl]-p^, such as fast Fourier transform, 
spectral methods and pseudospectral methods, or local ]T6|-p5t, such as finite difference, finite volume and finite 
element methods. Global methods are highly accurate but are cumbersome to implement in complex geometries 
and non-conventional boundary conditions. For example, a global method may converge slowly in a waveguide 
mode analysis due to irregular boundary conditions. In contrast, local methods are easy to implement for complex 
geometries and discontinuous boundary conditions. However, the accuracy of local methods is usually very low. 
There exists many problems in CEM which require both high computational accuracy and high numerical flexibility 
in handling complex geometries. These problems are characterized by a geometry which has a large domain size, i.e., 
the dimensions of the scatterer greatly exceed the wavelength of the incident wave. A typical example is the radar 
cross-section analysis of an entire airplane with an incident electromagnetic wave having a frequency of the order of 
ten GHz. To deal with such problems, it is desired to have a computational method that has both global methods' 
accuracy and local methods' flexibility. 

More recently, discrete singular convolution (DSC) algorithm was proposed as a potential approach for the computer 
realization of singular convolutions |26| , ^7[ |. Sequences of approximations to the singular kernels of Hilbert type, Abel 
type and delta type were constructed. Applications were discussed to analytical signal processing. Radon transform 
and surface interpolation. The mathematical foundation of the DSC algorithm is the theory of distributions ||2^ and 
the theory of wavelets. Numerical solutions to differential equations are formulated via the singular kernels of the 
delta type. By appropriately selecting parameters of a DSC kernel, the DSC approach exhibits controllable accuracy 
for integration and shows excellent flexibility in handling complex geometries and boundary conditions. Many DSC 
kernels, such as (regularized) Shannon's delta kernel, (regularized) Dirichlet kernel, (regularized) Lagrange kernel and 
(regularized) de la Vallee Poussin kernel, have been constructed p^ . Practical applications were examined for the 
numerical solution of the Fokker-Planck equation [p6|j2^ ] and for the Schrodinger equation . Another development 
in the application of the DSC algorithm is its use in computing numerical solutions of the Navier-Stokes equation and 
in structural analysis ||3^. In the context of image processing, DSC kernels were used to facilitate a new anisotropic 
diffusion operator for image restoration from noise Most recently, the DSC algorithm was used to integrate the 
(nonlinear) sine-Gordon equation with the initial values close to a homoclinic manifold singularity [ ^2| , for which 
conventional local methods encounter great difficulties and result in numerically induced chaos . 

The purpose of this paper is to study the computational philosophy of the DSC algorithm and to introduce the 
algorithm for computational electromagnetics (CEM). For the purpose of numerical computation, both bandlimited 
reproducing kernels and approximate reproducing kernels are discussed as sequences of approximations to the universal 
reproducing kernel, the delta distribution. A systematic treatment is proposed for handling a general class of boundary 
conditions. We explore the unified feature of the DSC algorithm for numerical approximation of differential equations. 
It is found that several conventional computational methods, such as methods of global, local, Galerkin, collocation, 
and finite difference can be derived from a single starting point. In particular, a Galerkin-induced collocation algorithm 
is discussed. A set of generalized finite difference schemes are shown to exhibit global-finite difference features at 
certain limit of DSC parameters. The present algorithm is shown to have controllable accuracy. The potential of 
the DSC algorithm for computational electromagnetics is explored by using three classes of problems, the eigenmode 
analysis of waveguide, the potential function analysis of electrostatics and the propagation of electromagnetic waves. 

This paper is organized as follows: In Section II, we study the computational philosophy of the DSC algorithm. 
A number of new DSC kernels are constructed as approximations to the universal reproducing kernel — the delta 
distribution. Approximation of functions and their derivatives is discussed. A systematic treatment of boundary 
conditions is proposed for being used in implicit schemes. The capability of the DSC algorithm is analyzed for solving 
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differential equations in Section III. The unified feature of the DSC algorithm is explored in the framework of the 
method of weighted residuals. The application of the DSC algorithm to CEM is introduced in Section IV. The 
utility and robustness of the proposed method is illustrated by a few numerical experiments. This paper ends with a 
conclusion. 



II. THE DISCRETE SINGULAR CONVOLUTION 



A. Approximation of singular convolution 

Singular convolutions (SC) are a special class of mathematical transformations, which appear in many science and 
engineering problems, such as Hilbert transform, Abel transform and Radon transform. It is most convenient to 
discuss the singular convolution in the context of the theory of distributions. The latter has a significant impact in 
mathematical analysis. Not only it provides a rigorous justification for a number of informal manipulations in physical 
and engineering sciences, but also it opens a new area of mathematics, which in turn gives impetus to many other 
mathematical disciplines, such as operator calculus, differential equations, functional analysis, harmonic analysis and 
transformation theory. In fact, the theory of wavelets and frames, a new mathematical branch developed in recent 
years, can also find its root in the theory of distributions. 

Let T be a distribution and ri(t) be an element of the space of test functions. A singular convolution is defined as 

/oo 
T{t~x)r]{x)dx. (1) 
'OO 

Here T{t — x) is a, singular kernel. Depending on the form of the kernel T, the singular convolution is the central 
issue for a wide range of science and engineering problems. For example, singular kernels of the Hilbert type have a 
general form of 

T(x) = ^, (n = l,2,...). (2) 

Here, kernel T{x) = ^ commonly occurs in electrodynamics, theory of linear response, signal processing, theory of 
analytic functions, and the Hilbert transform. When n — 2, T{x) = jj- is the kernel used in tomography. Another 
interesting example is singular kernels of the Abel type 

T{x) = ^, (0 < /3 < 1). (3) 

These kernels can be recognized as the special cases of the singular integral equations of Volterra type of the first kind. 
Singular kernels of the Abel type have applications in the area of holography and interferometry with phase objects 
(of practical importance in aerodynamics, heat and mass transfer, and plasma diagnostics). They are intimately 
connected with the Radon transform, for example, in determining the refractive index from the knowledge of a 
holographic interferogram. The other important example is singular kernels of the delta type 

r(x)-5(")(x), (n = 0,1,2,...). (4) 

Here, kernel T{x) = 5{x) is of particular importance for interpolation of surfaces and curves (including atomic, 
molecular and biological potential energy surfaces, engineering surfaces and a variety of image processing and pattern 
recognition problems involving low-pass filters). Higher-order kernels, T{^x) = ^("^(x), (n = 1, 2, . . .) are essential for 
numerically solving partial differential equations and for image processing, noise estimation, etc. However, since these 
kernels are singular, they cannot be directly digitized in computers. Hence, the singular convolution, ([I]), is of little 
numerical merit. To avoid the difficulty of using singular expressions directly in computer, we construct sequences of 
approximations (Tq,) to the distribution T 

lim Ta{x) — > T{x), (5) 

a — >ao 

where ao is a generalized limit. Obviously, in the case of T{x) = S{x)^ each element in the sequence, Ta{x), is a delta 
sequence kernel. Note that one retains the delta distribution at the limit of a delta sequence kernel. Computationally, 
the Fourier transform of the delta distribution is unity. Hence, it is a universal reproducing kernel for numerical 
computations and an all pass filter for image and signal processing. Therefore, the delta distribution can be used as a 
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starting point for the construction of either band-hmited reproducing kernels or approximate reproducing kernels. By 
the Heisenberg uncertainty principle, exact reproducing kernels have bad localization in the time (spatial) domain, 
whereas, approximate reproducing kernels can be localized in both time and frequency representations. Furthermore, 
with a sufficiently smooth approximation, it makes sense to consider a discrete singular convolution (DSC) 

F^{t)=Y,To.it-Xk)fixk), (6) 
fc 

where Fa{t) is an approximation to F{t) and {xk} is an appropriate set of discrete points on which the DSC (^) is 
well defined. Note that, the original test function ri{x) has been replaced by f{x). The mathematical property or 
requirement of f{x) is determined by the approximate kernel Tq,. In general, the convolution is required being Lebesgue 
integrable. In the rest of this paper, the emphasis is on the singular kernels of the delta type, their approximation, 
and numerical implementation. 



B. Singular kernels of delta type 

The delta distribution or so called Dirac delta function (6) is a generalized function which is integrable inside a 
particular interval but in itself need not have a value. Heaviside introduced both the unit step Heaviside function 
and the Dirac delta function as its derivative and referred to the latter as the unit impulse. Dirac, for the first time, 
explicitly discussed the properties of S in his classic text on quantum mechanics; for this reason S is often called Dirac 
delta function. However, delta distribution has a history which antedates both Heaviside and Dirac. It appeared in 
explicit form as early as 1822, in Fourier's Theorie Analytique de la Chaleur. The work of Heaviside, and subsequently 
of Dirac, in the systematic but informal exploitation of the step and delta function has made delta distribution familiar 
to physicists and engineers before Sobolev, Schwartz |2^ , Korevaar and others put it into a rigorous mathematical 
form. The Dirac delta function is the most important special case of distributions or generalized functions. There are 
three parallel descriptions for the theory of distributions. One description of distributions is to characterize them as 
an equivalence class, or as a generalized limit of various Cauchy sequences (fundamental sequences) and fundamental 
families as rigorously defined by Korevaar |]3^ . This approach is particularly convenient for the delta distribution. 
Another description is to formulate them as continuous linear functionals on the space of test functions as introduced 
by Schwartz The vector space of test functions is obtained from a class of test functions with compatible 

convergence or topology. The other description is based on generalized derivatives of integrable functions. Generalized 
derivatives are distributions rather than well-behaved functions. The first description is intuitive and convenient for 
various applications. The second description is particularly elegant and concise. It is also very convenient for higher 
dimensional applications. The third description is useful for certain practical applications involving derivatives and 
antiderivatives. These three descriptions are formally equivalent and are commonly used for describing not only for 
the delta distribution, but also for distributions in general. The use of many delta sequences as probability density 
estimators was discussed by Walter and Blum |3|] and others [^-^ . 

Definition 1. The delta distribution, or so called Dirac delta function is given as a continuous linear functional on 
the space of test functions, I?(— cx), cx)), 

/oo 
5(t> = 0(0). (7) 
-oo 

A delta sequence kernel, {6a{x)}, is a sequence of kernel functions on {—oo, oo) which is integrable over every compact 
domain and their inner product with every test function (j) converges to the delta distribution 



lim / Sa4> =< S,(j}>, (8) 

where the (real or complex ) parameter a approaches ao which can either be oo or a limit value, depending on the 
situation (such a convention for ao is used thorough out this paper). If ao represents a limit value, the corresponding 
delta sequence kernel is a fundamental family. Depending on the explicit form of Sa, the condition on (j) can be relaxed. 
For example, if Sa is given as 

a for < a; < l/a Q! = 1,2,--- , , 

otherwise ' ^ ' 

then Eq. (p|) makes sense for every (p in C(— oo, oo). 
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There are many delta sequence kernels arising in the theory of partial differential equations, Fourier transforms 
and signal analysis, with completely different mathematical properties. It is useful to have a classification of various 
delta sequence kernels for discussion, application and for new construction. Delta sequence kernels of positive and 
Dirichlet type have very distinct mathematical properties and can serve as the basis of a good classification scheme. 
In particular, there is a close relation between the delta sequence kernel of positive type and statistical distribution 
functions. In fact, every statistical distribution function can be regarded as an element of a delta sequence kernel of 
the positive type. An ordinary element of delta sequence kernel of Dirichlet type has the well-known feature of "small 
wave" . In other words it is readily related to the wavelet scaling function. Moreover, classifying delta sequence kernels 
according to Schwartz class or non-Schwartz class is also very useful for various physical and engineering applications. 
In particular, all physically realizable states, either in the sense of quantum mechanics or classical mechanics, belong 
to the Schwartz class Moreover, for the purpose of numerical applications to ill-posed problems, delta sequences 
of the Schwartz class are applicable to a wide class of functions and distributions. In the following two subsections. 
Delta sequence kernels of positive type and Dirichlet type are studied. 



C. Delta sequence kernels of positive type 

Definition 2. Let {Sa} be a sequence of kernel functions on (—00,00) which are integrable over every bounded 
interval. Wc call {da} a delta sequence kernel of positive type if 

1. J^^ Sa ^ 1 as a — > ao for some finite constant a. 

2. For every constant 7 > 0, + (J^ — > as a qq. 

3. Sa{x) > for all x and a. 

Example 1. Delta sequence kernel of impulse functions 

To approximate idealized physical concepts such as the force density of a unit force at the origin a; = 0, or a unit 
impulse at time a; = 0, a sequence of functions given by 

^ . X f a for < a; < 1/a a = 1, 2, • • • 

Oa\x) = i n +1, • (10) 

^ ' [0 otherwise ^ ' 

is a DSC delta sequence kernel provided a ^ 00. This is a commonly used density estimator in science and engineering. 
Example 2. Gauss' delta sequence kernel 
In the study of heat equation. Gauss' delta sequence kernel 

6a{x) ^ —^6-'='/^°'' for a^O (11) 
V27rck 

arises naturally as a distribution solution or so called weak solution. Gauss' delta sequence kernel has various 
interesting properties with regard to differentiability, boundedness and Fourier transforms, and it is used to generate 
the "Mexican hat" wavelet. 

Example 3. Lorentz's delta sequence kernel 

Lorentz's delta sequence kernel 

5a{x) ^ - " „ for a^O (12) 

is known for its role in representing the solution of Laplace equation in the upper half plane. It is commonly seen in 
integral equations involving the Green's function of the kinetic energy operator (in the momentum representation). 
It is also the expression for the line shape of various spectroscopies when the relaxation is an exponential one in the 
time domain. A generalized expression can be written as, 

1 a"a;"~^ 

5a,n{x) = 5— — 7T- for a ^ 0, and n > 1. (13) 

This includes Eq. ([l^ ) as a special case. 
Example 4. Landau's delta sequence kernel 
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In the discussion of convergence properties of polynomial approximations, Landau introduced a delta sequence 
kernel 



Ln{x) 



I ■ becomes a delta sequence kernel 



(a2 - x2) 



.2\n 



for n = 0, 1, 2, • • • and a > 0. 



W = { J" 



(x) for I a; |< a 
otherwise 



(14) 



(15) 



as n ^ 00. This is called Landau's delta sequence kernel. Wavelets generated from Landau's delta sequence kernel 

can be very useful for a sufficiently large n. 

Example 5. Poisson's delta sequence kernel family 

The function given by the summation of an infinite series 



1 

- + a cos(a;) 



cos{2x) 



27r(l - 2acos(x) +a2)^ 



where < a < 1 and (— oo < a; < oo), is called the Poisson kernel, which plays an important role in Poisson's integral 
formulae. Poisson's delta sequence kernel family is given by 



5a{x) 



Pa{x) for I a; |< TT 
otherwise 



(16) 



as a ^ 1. The Poisson's delta kernel family has a connection with the solution of Laplace equation in the unit disc 
(i.e., Dirichlct problem for the unit disc). 

Example 6. Fejer's delta sequence kernel 

The partial sum of the discrete Fourier series 



Dk{x) = - 

TT 



1 

sin[(fc 



cos(x) 
1 ^™l 



277 sin \x 



- cos(2x) H h cos(A:a;) 

fc = 0,1,2,- •• 



(17) 



is called a Dirichlct kernel. To improve convergence for proving a trigonometric approximation theorem, Fejer intro- 
duced the following arithmetic mean 



Fk{x) = i [Do{x) 



Di{x) + ■ ■ ■ + Dk-^{x)\ 



{\kx) 



2iTksh^{\x) 



00 < a; < 00. 



Then Fejer's delta sequence kernel is given by 

5a{x) -- 



Fa{x) for I .T |< TT for a = 0, 1, 2, • 

otlicTwise 



(18) 



(19) 



as a — * 00. Fejer's delta kernel has an important application in the theory of reproducing kernels. It also describes 
intensity pattern of light from a regular series of pinholes in optical physics. 

Example 7. Generalized Fejer's delta sequence kernel 

It is noted that Fejer's method of generating delta sequence kernel is very general. Essentially, a family of arithmetic 
means of delta sequence kernels is still a delta sequence kernel. The resulting delta sequence kernel can be called a 

delta sequence kernel of Fejer type. For instant, in a similar treatment using Dirichlet's continuous delta sequence 
kernels (see next subsection), one obtains the following Fejer's continuous delta sequence kernel 



- , , 2 sin^(aa;) ^ „ 
5a{x) = Va; e R. 



ax^ 



(20) 
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Obviously, this is well-defined on the real line. This expression is related to the intensity of light diffracted by a 
uniform slit. 

Example 8. Delta sequence kernels generated by dilation 

Let p G L^{R) be a non-negative function with J p{x)dx ~ 1, dilation of p is given by 

p^{x)^-p{-) (a>0) (21) 
a a 

leads to a delta sequence kernel, pa ^ S, as a ^ 0. 

Physically, p can be regarded as a statistical distribution function. This is a general procedure and Examples 2 and 
3 fit into this structure. Examples 1 and 6 can be expressed in this form by appropriate modifications (by replacing 
a with (3 = and then letting /3 ^ 0). 



D. Delta sequence kernels of Dirichlet type 



Definition 3. Let {5a} be a sequence of functions on (—00,00) which are integrable over every bounded interval. 
We call {5a} a delta sequence kernel of the Dirichlet type if 

1. J^^ (5q — > 1 as a — > ao ^oi some finite constant a. 

2. For every constant 7 > 0; ^ -I- (5a as a ao, 

3. There are positive constants Ci and C2 such that 

|<5.(x)|<^-fC2 

for all X and a. 
Example 1. Dirichlet kernel 

The most important example of a delta sequence kernel of Dirichlet type is Dirichlet kernel 

c ( \ _ \ Da{x) for|a;|<7r for a^O, 1,2,--- , . 

Oo.{x) -<y^ otherwise ^ > 

where Da is the Dirichlet kernel given by Eq. (p7|). Dirichlet's delta sequence kernel plays an important role in 
approximation theory and is the key element in trigonometric polynomial approximations. In fact, it is an exact 
reproducing kernel for bandlimited, periodic, functions. Physically, it describes the diffraction of light passing 
through a regular series of pinholes in which the fcth pinhole's contribution is proportional to e**^. 
Example 2. Modified Dirichlet kernel 

Sometimes there is some advantage in taking the last term in Da with a factor of i; 

Da{x) ^ Da~ cos{ax) 
Zn 

- a = 0,1,2,.... (23) 



2tt tan(ix) ' 

This is the so-called modified Dirichlet kernel. The difference Da — D*^ tends uniformly to zero on (— tt, tt) as a ^ 00. 
They are equivalent with respect to convergence. 
The expression given by 

Arr.^-/^"(^) foi'|a;|<7r for a = 0, 1, 2, . . . 

Oo.{^) - I otherwise ^ > 

is a delta sequence kernel of Dirichlet type as a ^ 00. 
Example 3. Lagrange kernel 
Lagrange interpolation formula 
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i=k+M ^ ^ 



LmA^)^ n ^^-3^' (^^^1) (25) 

i=k-M,i^k 

is defined on an interval (a, h) with a set of 2M + 1 ordered discrete points, 

{^i}'i=k-M ■ ^k-M = a < Xfe-M+i < ■ • • < a;fe < • • • < Xfc+Ai = b. (26) 
It converges to tlie delta distribution as 

a — > —00, b 00 and sup | Xi — a;j |— > 0. (27) 

Obviously, these limits imply M 00. Since the delta distribution has only a point support, the Lagrange interpo- 
lation formula is a delta sequence 

c / ^ { LM.k{x) for a < X < & for A'f = l,2,--- , - 

= otherwise (^8) 



as, M 00 (to qualify as a delta sequence of the Dirichlet type, M > 2 is required in Eq. ( |2g 
Example 4- Interpolative delta sequence kernel 

Let {Sn} be a sequence of functions converging to the delta distribution and let {x^Iq be n + 1 zeroes of a Jacobi 
polynomial in (a, 5). 

\^ y) iii=oyy ^'■i 

is a delta sequence kernel as n 00. This follows from the fact that J An{x,y)f{y)dy are approximations to the 
Lagrange interpolation formula. 

Example 5. The de la Vallee Poussin delta sequence kernel 

The de la Vallee Poussin kernel is given by 



1 " 

P I ^ ^ 



k—n—p 



1 1 



Ecos kx A — > 
TT ^ 

fc=l fc=l 

sin[(2n+ 1 -p)|]sin[(p+ l)f] 
27r(p+l)sin2(f) 



fc 



cos[(n — p + fc)a;] 



p+l_ 

, p = Q,---,n; n = 0, 1,---, (30) 



where Dk{x) are the Dirichlet kernels given by Eq. (0). It is interesting to note that de la Vallee Poussin kernel 
reduces to the positively defined Ferer's kernel Fn+i{x) when p — n. The de la Vallee Poussin delta sequence kernel 
is given by 

, , J F„_p(x) for I a; I < TT for p = 0, • • • , n; n = 0, 1, • • • , , . 

"n.pia;; - I g otherwise ^'^ ' 

as n,p ^ 00. The de la Vallee Poussin delta sequence kernel is of Dirichlet type when p < n. 
A simplified de la Vallee Poussin kernel given by 

^ , , 1 cos(Q!a:) — cos(2aa;) , , 

na x'^ 

is found to be very useful numerically [p6| . 

Example 6. DSC kernels constructed by orthogonal basis expansions 
Let {tpi} be a complete orthonormal L'^{a,b) basis. Then 

n 

6n{x,y) ^^'^p^{x)■^p^{y), x,y e {a,b) (33) 

i=0 
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are DSC delta sequence kernels. In the case of trigonometric functions, we again obtain the Dirichlet kernels given in 
the Examples 1 and 3. Hermite function expansion is given by 

(5„(a;)=exp(-x2)^ ( — ) ^^H2k{x), e R (34) 



where H2k (x) is the usual Hermite polynomial. Note that the Hermite's kernel in Eq. (^J) has a different form from 
Eq. (p^). This is because it is evaluated at a; = in the series expansion. 
Example 7. Shannon's delta sequence 

Shannon's delta kernel (or Dirichlet's continuous delta kernel) is given by the following (inverse) Fourier transform 
of the characteristic function, x\-a q1: 



1 f°° 
sin(a.x) 



TTX 

Alternatively, Shannon's delta kernel can be given as an integration 



(35) 



1 

Sa{x) = - cos{xy)dy, (36) 
Ti" Jo 

or as the limit of a continuous product 

r / N T Q; TT f \ T 1 sm[ax) 

SJx) = hm - cos{-rx) = hm — rj . ,\ ' . (37 

^ ' w^oo TT J-J^ ^2*^ ' JV^oo 2^7rsin(^a;) ^ ' 

Numerically, Shannon's delta kernel is the most important, because of its property of being an element of the 
Paley- Wiener reproducing kernel Hilbert space B^. 

m = r fiy f^f'^'f dy^ v/ e bi, (38) 

J-oc T^lx - y) 

where V/ G B^. indicates that, in its Fourier representation, the function / vanishes outside the interval [— 7r,7r]. 
The Paley- Wiener reproducing kernel Hilbert space B'^ is a subspace of the Hilbert space L^{R). It is noted that the 
reproducing kernel Hilbert space is a special class of Hilbert space. For instance, the space (i?) is not a reproducing 
kernel Hilbert space. 

Example 8. Generalized Lagrange delta kernel 

Shannon's delta kernel can be derived from the generalized Lagrange interpolation formula 

where G{x) is an entire function given by 

G[x) = {x- xo) n (l - ^) - ^) ' (40) 

and G' denotes the derivative of G. For a function bandlimited to B, the generalized Lagrange interpolating formula 
Sk [x) of Eq. ( |3^ ) can provide an exact result 

fix) = J2fiyk)Sk{x), (41) 

kez 

whenever the set of non-uniform sampling points satisfy 



sup 

kez 



Xk-^ 
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where the symbol Z denotes the set of all integers. This is called the Paley and Wiener sampling theorem in the 
literature. 

If {xk}kez are limited to a set of points on a uniform infinite grid {xk — kA — —x-k), Eq. ( ^0| ) can be simplified 

oo 

G(-)=- n (i-^) (43) 

fc— — oo, fc^O 

(45) 



X 

k=l ^ 

. sin 



n 



Since G'{xk) reduces to 



on a uniform grid, Eq. (BS) gives rise to 



G'{xk) = {-If (46) 



G{x) _ (-I)'' sin fa; 
G'{xk){x-Xk) ~ f{x-kA) 



sinf (a; - Xk) 

^{X-Xk) 

) ,v „ 



(48) 



Obviously, '^™a(^ ^fc) jg approximation to the delta distribution 



sinf (x-Xfc) 

l^^o f -^(^-^'c). (49) 

In fact, the generalized Lagrange interpolation formula directly gives rise to the delta distribution under an appropriate 
limit 

lim S'fc(x) = lim ^}f^ -~^&{x-Xk), (50) 

maxAx-»0 m&y. l^x^Q G {Xk)\X — Xk) 

where max Aa; is the largest Ax on the grid. 

E. Connection to wavelets 

The DSC approximation to the delta distribution is closely related to the theory of wavelets and frames. Mathe- 
matically, wavelets are functions generated from a single function by applying dilation and translation. They form 
building blocks for some spaces, such as L^{R), whether as a frame or as an orthonormal basis. Such building blocks 
are computationally important when they have certain regularity and localization in both time and frequency do- 
mains. Physically, the wavelet transform is a mathematical technique that can be used to split a signal into different 
frequency bands or components so that each component can be studied with a resolution matched to its scale, thus 
providing excellent frequency and spatial resolution, and achieving computational efficiency. 

Shannon's wavelet is one of the most important examples and its scaling function is the Shannon's delta kernel, 

a; = . (51) 

TTX 

As a delta kernel, it is normalized 

,^(0) = J (j){x)dx = 1, (52) 
and its Fourier transform is given by the characteristic function (t>{Lu) = X[- 1/2.1/2)- It is easy to see that 
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+ n) = l (53) 

71 — — OO 

and 

oo 

^ \${iv + n)\'=l. (54) 



n— — OQ 



Equation ( |5^ ) is a consequence of orthonormality. In fact, the sequence of functions {(/)(a; — 7i)}^_oo orthonormal. 
Shannon's mother wavelet can be constructed from the Shannon's delta kernel (Shannon's wavelet scaling function) 

sin 2t:x — sin nx 

V[x) = , (55) 

TTX 

with its Fourier expression 

i'i^) = X[-i,i] i'^) - Xhi/2a/2] ('^)- (56) 
This is recognized as the ideal band pass filter and it satisfies the orthonormality conditions 

oo 

V>(w + n) = 1 (57) 

n— — OO 

and 

oo 

Y \kc^ + n)f=l. (58) 



n— ~oo 



Technically, it can be shown that a system of orthogonal wavelets is generated from a single function, a "mother" 
wavelet ip, by standard operations of translation and dilation 

?/;„„(x) = , m,neZ. (59) 

A family of Shannon's wavelet scaling functions {'>jjmn{x)}n.mez span a series of orthogonal wavelet subspaces 
{Wm}m.GZ satisfying 

^Wm^L^R). (60) 

mez 

Alternatively, a family of Shannon's wavelet scaling functions {4>mn{x)}n,mez arc constructed from a single Shannon's 
delta kernel 

</>mn(a;) = 2"^(/) - , m,neZ. (61) 

They span a series of nested wavelet subspaces {Vm\mez, each corresponds to a different resolution 

••• C VLi C C C ••• C i^(i?). (62) 

This nested structure provides the conceptual basis for the wavelet multiresolution analysis. 

From the point of view of signal processing, the Shannon's delta kernel 0q corresponds to a family of ideal low pass 
filters, each with a different bandwidth 



<Pa{x) = . (63) 



Their corresponding wavelet expressions. 



, , sin 2a'KX — sin airx 

{x) - , 64 

TTX 

are band pass filters. However, the Shannon's wavelet system is seldom used in real applications because it requires 
infinitely many data points. In the next subsection, we discuss a practical approach for generating powerful filters 
from Shannon's delta kernel. 
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F. Regularization 



Both (l){x) and its associated wavelet play a crucial role in information theory and the theory of signal processing. 
However their usefulness is limited by the fact that (j){x) and ip{x) are infinite impulse response (IIR) filters and their 
Fourier transforms 4>{i^) and tp{to) are not differentiable. From the computational point of view, (j){x) and tp{x) do 
not have finite moments in the coordinate space; in other words, they are de-localized. This non-local feature in 
the coordinate is related to its bandlimited character in the Fourier representation by the Heisenberg uncertainty 
principle. 

According to the theory of distributions, the smoothness, regularity and localization of a temper distribution can 
be improved by a function of the Schwartz class. We apply this principle to regularize singular convolution kernels 

^aix) = Ra{x)(l){x), (a > 0), (65) 

where Ra- is a regularizer which has properties 

lim R^{x) = 1 (66) 

a — ^oo 

and 



Ra{0) = 1. (67) 



Here Eq. (|66|) is a general condition that a regularizer must satisfy, while Eq. (|6^) is specifically for a delta regularizer, 
which is used in regularizing a delta kernel. Various delta regularizers can be used for numerical computations. A 
good example is the Gaussian 



Ra{x) = exp 



x^ 
'2^2 



(68) 



Gaussian regularizer is a Schwartz class function and has excellent numerical performance. However, we noted that 
in certain eigenvalue problems, no regularization is required if the potential is smooth and bounded from below (e.g., 
the harmonic oscillator potential \x'^). 

Immediate benefit from the regularized Shannon's kernel function Eq. ( |65| ) is that its Fourier transform is infinitely 
differentiable because the Gaussian is an element of Schwartz class functions. Qualitatively, all kernels of the Dirichlet 
type oscillate in the coordinate representation. Shannon's kernel has a long tail which is proportional to ^, whereas, 
the regularized kernels decay exponentially fast, especially when the a is very small. In the Fourier representation. 
Shannon's kernel is an ideal low pass filter, which is discontinuous at w = i. In contrast, all regularized Shannon's 
kernels have an "optimal" shape in their frequency responses. Of course, they all reduce to Shannon's low pass filter 
at the limit 

sin TTX sin TTX , , 

lim <Po-(x) = lim e ^ = . (69) 

a-^oo cr— >oo TTX TTX 

Quantitatively, one can examine the normalization of ^a{x) 



$a(0) = / ^a{x)dx 

(-1)^= /7r(T\"^" 



2vra^ 



fc=0 



fc!(2fc + l) VV2/ 



(70) 



By means of the error function, erf(z) = e * dt, Eq. ( |70[) can be rewritten as 

*,(0) = erf(^| 



oo 



e 2 ; e dt 

'^(^ Jo 



1 _ erfc I ^ ) , (71) 
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where erfc(z) is the complementary error function. Note that for a given cr > 0, erfc(^) is positive definite. Thus, 

4>(j(0) is always less than unity except at the hniit of cr —> oo. 

In fact, ^cr{x) does not really satisfy the requirement, as given by Eq. (^^, for a wavelet scaling function. However, 
when we choose a ^ V^/tt, which is the case in many practical applications, the residue term, erfc(-^), approaches 

zero very quickly. As a result, $ct(0) is extremely close to unity. Therefore, we call the regularized Shannon's delta 
kernel $0- a quasi-wavelet scaling function. 



G. Discretization 



For the purpose of digital computations, it is necessary to discretize various delta kernels. To this end, we should 
examine a sampling basis given by the Shannon's delta kernel 

Sk[x) = K{x,Xk) = Vfc e Z. (72) 

n[x - Xk) 

This sampling basis is an element of the Paley- Wiener reproducing kernel Hilbert space. Hence, it provides a discrete 
representation of every (continuous) function in that is 

/(x) = ^/(a;fe)5,(a;), V/ G B^. (73) 

kGZ 

This is recognized as Shannon's sampling theorem and it means that one can recover a continuous bandlimited 
function from a set of discrete values. Equation (|7^) is particularly important to information theory and the theory 
of sampling because it satisfies the interpolation condition 

where S^.k is the Kronecker delta function. Note that Shannon's delta kernel is obviously interpolative on Z. Com- 
putationally, being interpolative is desirable for numerical accuracy and simplicity. 
On a grid of arbitrary spacing A, Shannon's sampling theorem can be modified as 



/(^)-E/(^^) ti^rjw V/eM. (75) 



This suggests that we can discretize the regularized Shannon's delta kernel as 

^ / ^ smf{x-Xk) .ii^fili 

^{X-Xk) 

It is noted that if A is chosen as the spatial mesh size (this is, in general, not required in signal and image processing), 
^a;A{x — Xk) retains the interpolation property, 

^a.Ai^rn ~ Xk) = Sm.k- (77) 

This is of particular merit for numerical computations. 

In practical applications, Eq. (|7^) can never be realized because it requires infinitely many sampling points. There- 
fore, it is both necessary and convenient to truncate the infinite summation in Eq. ( [75|) to a finite (2M + 1) summation 

M 

/(x)« ^-A{^-^k)f{xk), (78) 

where S^^^ix — Xk) is a collective symbol for any (regularized) delta kernel. The truncation error is dramatically 
reduced by the introduction of a delta regularizer. A rigorous proof of this has been given by Qian and Wei pC| ]. 

The discretization the of the Dirichlet kernel is not as straightforward as is for the Shannon's kernel. However, it 
can be carried out according to the following Dirichlet sampling theorem: 
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Theorem // an function f{x) satisfies the Dirichlet boundary condition and is periodic in T and bandlimited 
to the highest (radial) frequency 2tiL/T, it can he exactly reconstructed from a finite set of 2L + 1 discrete sampling 
points 



L 

E 

fe=-L 



f{xk) 



(2L + l)sin 



A 2L + 1 



(79) 



where A — T/(2L + 1) is the sampling interval and Xk = fcA are the sampling points. 

Note that the kernel in Eq. ( [79|) differs form that in Eq. This follows from a change in the variable 

X —> -^jT+iy' ^^'^ / '^•^ ~^ X] 'Sy 2L+1 '^'^^ Dirichlet kernel is a reproducing kernel for bandlimited periodic 
functions. Therefore, Eq. (|7^) should be the most efficient kernel for numerical computations under the aforementioned 
conditions. However, to facilitate the Dirichlet kernel in an unbounded computational domain, we use the following 
regularized discrete expression for the Dirichlet kernel: 



sm 



[{l + \){x-x')\ sin[i(x~Xfc)] 



27r sin[^(a; — x')] 



(2L + l)sin 



A 2L+1 



■exp 



Xk)'^ 



(80) 



Like the regularized Shannon's kernel filter, the present regularized Dirichlet kernel filter has rapid decay properties. 
In comparison to Shannon's kernel, the Dirichlet kernel has one more parameter L which can be optimized to achieve 
better results in computations. Usually, we set a sufficiently large L for various numerical applications. A regularized 
discrete expression for the modified Dirichlet kernel is 



sm 



[[I + ^] (x - x')] sin[f{x-Xk)] 



27rtan[i(a; 



X')] 



(2L + l)tan 



7r x — Xk 
A 2L + 1 



■exp 



(X - Xfc)^ 

2ct2 



(81) 



Obviously, the regularized Dirichlet kernel reduces to the regularized Shannon's delta kernel when L is sufficiently 
large 



lim ■ 

L — >oo 



lim 

L — ^oo 



sin [^{x - Xk) 



TT_ X~Xk 

A 2L + 1 



[2L+ l)sin 

sin [^{x - Xk)] 



■exp 



{x - XkY 



(2L+ l)tan 



sin^(a; - Xk) 



TT_ X~Xk 

A 2L + 1 



■exp 



2ct2 



{x - XkY 



2a2 



Xk) 



(82) 



The discretization of the de la Vallee Poussin kernel is given by 

1 cos[a(x — x' )] — cos[2a(a:: — x' )] 



[x — x' )2 
2 cos -^(a; — Xk) — cos "^{x 



Xk) 



[fix-Xk)]' 



-exp 



{x - Xk)'^ 

2ct2 



(83) 



where A = | A. Since vr/A is proportional to the highest frequency which can be reached in the Fourier representation, 
the A should be very small for a given problem involving very oscillatory functions or very high frequency components. 
It is noted that by definition, the Lagrange interpolation formula 

k+M 



i=k-M,i^k 



is already discretized. However, its regularized forms 

k+M 

Saix-Xk)= Yl 

i=k-M,i^k 



Xk Xi 



Xk 



exp 



(a: - XkY 

2(72 



(85) 



are very good low pass filters. Unlike the generalized Lagrange interpolation formula, both expressions ( |84| ) and ( |85| ) 
are compactly supported kernels. As low pass filters, they require only only a finite number of signals and their 
Fourier transforms have smoothened shoulders as those of regularized Shannon's delta kernels. 
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H. Approximation of derivatives 



For the solution of differential equations, an approximation to the derivatives is required. Such an approximation 
can be constructed by using DSC kernels of the delta type with n 7^ 0. Let us consider a one-dimensional, nth order 
DSC kernel of the delta type 



(x-Xk), (n = 0, 1,2, • • •). 



(86) 



Here Sl^]^{x — Xk) — 6cr.A{x — Xk) is the DSC delta kernel described in Eq. (|78|). The higher order derivative terms 
— Xk) are given by differentiation 



— ) 5,^a{x - Xk) 



(87) 



These derivatives can be regarded as high pass filters. The filters corresponding to the derivatives of Shannon's kernel 
decay slowly as x increases, whereas, the regularized filters are Schwartz class functions and have controlled residual 
amplitudes at large x values. In the Fourier representation, the derivatives of Shannon's kernel are discontinuous at 
certain points. In contrast, all the derivatives of regularized kernels are continuous and can be made very close to 
those of Shannon's, if desired. 

The differentiation in Eq. ( ^7|) can be analytically carried out for a given 5c^/\{x—Xk)- For example, if 6rj^/\(x — Xk) = 

(x-xj.)2 

, we have for x ^ Xk 



-(x-Xk) 



5':^A^X-Xk) 



COS ^ [x 


- Xk) 


{x- 


Xk) 


sin -J- {x 


- Xk) 




xky 


sin -J {x 


- Xk) 



exp 
exp 
exp 



(x-Xk)^ 

2(72 

jx- Xk)^ 
2(72 

{X - Xk)^ 
2cr2 



^sin^(x - X k) 

{X - Xk) 
COS^(x - Xk) 

[X - Xk)^ 
COSf (X- Xfc) 

2 ^^-^5 exp 



■ exp 



{X - Xk)' 



sinf{x-Xk) 
iix-Xk)^ 
sinf (a; - Xk) 



exp 



^(X - Xfc)t72 

sin 3- (a; - Xk) 



exp 



2a2 

{X - Xfc)2 
2(72 

jx-Xk)"^ 

2(72 

{x ~ Xk)^ 
2(72 

{X - Xk)'^ 
2(72 



i(74 



(a; - Xk) exp 



{X - Xfc)^ 

2a2 



(89) 



^cos^{x - Xk) ( {x-Xkf 
— exp ' 



{x - Xk) 

fsin J(x~ Xfc) 
3 , 7^ exp 

[X - Xk)^ 

f sinf(x- Xfc) 
3-^^^ ^ exp 



2(72 

{X - Xk)'^ 
2(72 

(X - Xk)^ 
2(72 



COS^{x-Xk) f {x~Xk)^ 

6 — — — exp ' 



(a; - Xk)^ 



2a2 
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and 



,C0S^{X - Xk) f {x-Xkf 



[X — Xkjcr \ Za'^ 

(x - Xk)cos^{x - Xk) f {x-XkY' 
+ 3 exp 



sin^ix - Xk) f {x-Xk f 



^ sin^jx-Xk) ( {x-Xkf 
\{x-XkYcj^^''^\ 2a2 

{x - xkf sill ^ (.r - Xk) ( [x- Xkf 



exp 0,2 h (90) 



- Xfc) = . - exp 



^sinf(a;-Xfc) / {x - Xkf 
+ -'^^ — — ^ exp ' 



(a; — Xk) \ 2o-2 

, ^cos^(a;-a:fc) / (^^^ 

+ a2 ^""Pl, 2a2 

lo jsinf (x-a:fc) ^ ^ / (a; - a;fe)2 
7 7^ exp I 



{x — Xk)^ V 20-2 

Jsinf (a;- Xfc) / (a; -2:^)2 
- 6-^^— — - — -^5 — exp ' 



(a; — Xk)cr'^ \ 2a'^ 

^{x - Xk)smf{x - Xk) ( 
- ^- exp - 



_ COsf(x-Xfc) ( (x-Xfc)2 

z4 . exp 



2a2 



(a; - Xk)^ ^ V 2a2 

_ .p cosf (x-Xfe) ( (x - a:fc)2 

(x-Xfe)2a2 ^""P^, 2a2 

(a; - a;fe)2 cos f (a; - a;fc) / {x - Xk)"^ 
— 4 exp — - — 

Cr6 V 20-2 

+ 24— TT^^ exp 



lix-Xk)^ "V 2a2 
^^^ sini(x-x.) ^ /_(x-x.)2 



f{x-Xk)^a^ '^\ 2a2 

sinf(a; - Xfc) / (x-Xfe)2 



A 

At X = Xfc, it is convenient to evaluate these derivatives separately 



f(.T-Xfe)a4 ^ " V 2^2 

^ (x-Xfc)sinf (x-Xfc) / (x - Xfc)2 
^ ^"^y 2a2 

^ (X-X.)3y(x-X.) ^^ ^g^^ 



<5i!i(0)=0 (92) 

^1(0) = (93) 
5(^(0) =0 (94) 
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and 



(0) 



1 15 + lOf^g^ + 

5 cr4 



(95) 



Similar expressions for other DSC kernels described in the last subsection can be easily derived. The performance 
of a few DSC kernels for fluid dynamic computations and structural analysis was given in Ref. Note that the 

differentiation matrix in Eq. ( |8^ ) is generally banded. This has a distinct advantage in large scale computations. 

For numerical computations, it turns out that the approximate reproducing kernel has much less truncation errors 
for interpolation and numerical differentiation. Qian and the present author | |40| ] have recently given the following 
theorem for truncation errors. 

Theorem Let f be a function f E L^{R) H C'^{R) and handlimited to B, {B < A is the grid spacing). For a fixed 
t E R and a > 0, denote g{x) — f{x)Hk{-^^^), where Hk{x) is the kth order Hermite polynomial. If g{x) satisfies 



g'{x) < g{x^ 



(x-t) 



(96) 



for x>t+ {Ml - 1)A, and 



g'{x) > g{x 



{x-t) 



(97) 



for X < t — M2A, where Mi, M2 G Af, then for any s G 

"sin^(i- nA) 



-(t-nA) 



exp(- 



{t - nAf 



is) 



27ra(i-B)exp(^^(^^) 



\f{t)\\LHR)Y. 



1 w ' -"1^ \ 



i+j+k=s i!fe!A'-i(V2(T)'=((Afi-l)A)J + i 



exp(^) 



-M2 A 
V2cr 



s i!fc!A'-i(V2(T)'=(M2A)3 + i 



exp( 



(M2A)2 
2a2 



(98) 



where superscript, (s), denotes the sth order derivative. 

The proof and detailed discussion (including a comparison with the truncation errors of Shannon's sampling theo- 
rem) are given in Ref. | |4C| ] and are beyond the scope of this paper. This theorem provides a guide to the choice of M , 
a and A. For example, in the case of interpolation (s = 0), if the L2 norm error is set to 10^'' (77 > 0), the following 
relations can be deduced from Eq. (Bq) 



and 



-{tt - BA) > y/4~61r], 



M , 

— > V4-61?7, 
r 



(99) 



(100) 



where r = cr/A (The choice of a is always proportional to A so that the width of the Gaussian envelope varies with 
the central frequency). The first inequality states that for a given grid size A, a large r is required for approximating 
high frequency component of an function. The second inequality indicates that if one chooses the ratio r = 3, 
then the half bandwidth M ^ 30 can be used to ensure the highest accuracy in a double precision computation 
(77 = 15). However, for lower accuracy requirement, a much smaller half bandwidth can be used. In general, the value 
of r is proportional to M . An appropriate value of M is determined by the accuracy requirement. This theoretical 
estimation is in excellent agreement with an earlier numerical test El[. 
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I. Treatment of boundary conditions 



While solving a set of differential equations, boundary conditions should be satisfied. In a global method, the kernel 
must be constructed in an adaptive manner near the boundary. The major drawback of such an algorithm is its 
difficulty in the construction of adaptive kernels for problems involving complex geometries and boundary conditions. 
Thus, global methods have been relatively less successful in dealing with these problems compared to their advantage 
in solving problems with simple geometries and boundary conditions. In contrast, the DSC algorithm utilizes a 
completely different philosophy for the kernel construction — the differentiation kernel is the same everywhere and is 
translation invariant on the grid. Therefore, the DSC algorithm is very flexible in dealing with a variety of boundary 
conditions and geometries. Moreover, the DSC treatment of boundary condition has its mathematical justification. 
In fact, for a continuous function, a derivative at a point exists if and only if both the left and the right derivatives 
exist and are equal. Therefore, with a finite computational domain, a boundary condition involving differentiation 
(such as the Neumann boundary condition) does not make sense. This is an ill-posed problem from the point of view 
of mathematical differentiation. The DSC approach to this problem is to extend the domain of definition for the 
system so that the "boundary condition" is well defined (i.e., differentiation is well-defined right on the boundary). 
The philosophy behind DSC is that, the original singular convolution has to be recovered from the DSC at the limit 
of Aa; —^ 0, everywhere in the computational domain, including boundaries. Therefore, at a boundary, a fictitious 
domain is required to ensure that the boundary condition is exactly satisfied at the continuous limit. 

In an explicit treatment, boundary conditions are easily implemented in the DSC algorithm by appropriate boundary 
extensions, which were discussed in Ref. In the present work, we consider implicit cases, where the boundary 

conditions are to be satisfied by a set of linear algebraic equations. 

We assume that a general boundary condition at the left boundary is given by 



N 



(101) 



n=0 



where Kn is a constant and xq is the boundary point. In the DSC treatment, we make an assumption for the relation 
between the inner nodes and the outer nodes on the left boundary 



f{x-i) - f{xo) = a, [{f{x,) - /(xo)] , 
where parameters a^, i = 1, . . . , M are to be determined. After a rearrangement, we obtain 

/(x_,) = aj{x,) + (1 - a,)f{xo), z = 1, 2, . . . , M. 



(102) 
(103) 



According to Eq. (p?]), we approximate the nth derivative and second derivative at the left boundary by 

M 



i=-M 

M 



(104) 



where C" are coefficients given by a DSC kernel, e.g., — <5^"^(xo — Xi). Therefore, the boundary condition (101) 
on the left boundary is given by 



N 



Kofixo) + E [CSfixo) 
n=l 

M 

i=l 
M 



1=1 



= 0. 



(105) 



In general, Eq. (10£) may not have an exact solution for implicit schemes. However, for a class of boundary conditions 
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TV 



/(xo) = 0, ^7^„/(")(xo) = 0, 



n=l 



the exact solution is given by 



(106) 



~,2;-i 



i = l,2,...,M. 



(107) 



The following a few special cases are important in practical applications. 
One case is the clamped edges, i.e., the boundary conditions require 

/(xo) = 0, fW{xo)^0. 

These are satisfied by choosing = 1, i = 1, 2, . . . , M. This is the so call ed s ymmetric extension 
Another case is simply supported edges, i.e. the boundary conditions (|101|) reduce to 



/(xo) = 0, /(2)(a;o) = 0. 

These are satisfied by choosing a.; = —1, i = 1,2,..., M. This is the so called anti- symmetric extension j26| . 
The other special case of Eq. (101) is called transversely supported edge and the boundary conditions are , 



/(a;o) = 0, f'^^Hxo) + Kif^'Hxo) ^ 0. 



In this case, Eq. (105) is reduced to 



M 



^[(1 + a^Cf +K,{1~ a.,)Cl]f{x,) = 0. 



i=l 



One way to satisfy Eq. (Ill) is to choose 



^^^'+^? ,-12 M 



(108) 

(109) 

!;iven by 
(110) 

(111) 
(112) 



Expressions for the right boundary can be derived in a similar way. 



III. UNIFIED FEATURES FOR SOLVING DIFFERENTIAL EQUATIONS 

To solve a differential equation, one can start, either by approximating the original differential operator or by 
approximating the actual solution of the differential equation while maintaining the original differential operator. 
The latter is accomplished by explicitly defining a functional form for approximations. Let us assume that the 
differential equation has the general form 

Cu{x)^f{x), xen, (113) 

where >C is a linear operator and u{x) is the unknown solution of interest. Here f{x) is a known source term, f2 
denotes the domain over which the differential equation applies. 

The approximate solution is sought from a finite set of N DSC trial functions of a given resolution a, denoted by 
with M being the half width of support of each element. Here cr is a regularization parameter for improving the 
regularity of the set. The case of regularization free is easily obtained by setting a 00. Elements of the set S^'^ 

can be exphcitly given by {4'a,a-ii 4>a^.<j-2^ 4>acr-N}- -^^^ ^ given computational domain, the resolution parameter a 
is determined by N. 

An important property of the DSC trial functions {(l>a(j-k} that when the trial function is free of regularization, 
each member of the set is a reproducing kernel at highest resolution 

lim <<,;fe,r7>=^(a;fe), (114) 
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where < •, • > denotes the standard inner product. In fact, if an appropriate basis is used for (f> and the Hmit on a is 
taken, of each resolution can be a reproducing kernel for functions bandlimited to appropriate sense. In general, 
we require the low pass filter property that for given a 7^ 0, ct 7^ and M ^ Q 



< 



(115) 



This converges uniformly when the resolution is refined, e.g., a — > 00. A few examples of such DSC trial functions 
are given in Refs. and and many more examples are constructed in the previuos section. Equations (114) 
and (115) are special requirements satisfied by the DSC kernels of delta type 

In the present DSC approach, an approximation to the function of interest u{x) can be expressed as a linear 
combination 



N 

E 

fc=i 



Ua,a;k4>a.(j\k{^)j 



(116) 



where x is an independent variable and Ua.a-.k is the desi red DSC approximation to the solution at point Xk- This 
structure is due to the DSC trial function property (11£) and it dramatically simplifies the solution procedure in 
practical computations. 

In this formulation, we choose the set S^'^^ a priori, and then determine the coefhcients {Ua,a;k} so that U^'^'^{x) 
is a good approximation to u{x) . T o determine Ua.a-.k, we minimize the amount by which U^'^^{x) fails to satisfy 
the original governing equation (113). A measure of this discrepancy can be defined as 



R 



N,M( 



ix)^CUZf{x)-f{x), 



(117) 



where R„' ^ (x ) is the residual for a particular choice of resolution, regularization and half width of the support. Note 



that Eq. (117) is similar to the usual statement in the method of weighted residuals. However, the appr oxim ation 
Ua'a^{x) is constructed by using the DSC trial functions, (p'^^-ki^)^ in the present treatment. Let Eq. (113) and 
its associated boundary conditions be well-posed, then there exists a unique solution u{x) which generally resides in 
an infinite-dimensional space. Since the DSC approximation U^'^^ is constructed from a finite-dimensional set, it is 
generally the case that U^;^^{x) ^ u{x) and therefore R^;^'' {x) ^ 0. 

Galerkin. We seek to optimize R^'^'^{x) by forcing it to zero in a weighted average sense over the domain Vl. A 
convenient starting point is the Galerkin 



R^f{x)K,\^,.^,{x)dx = 0, G SZ;:^ , 



M' 



-,N',M' 



(118) 



^N' ,M' 



where the weight set S^, can be simply chosen being identical to the DSC trial function set S'4 ^ . We refer to 



Eq. (lis) as a DSC-Galerkin statement. 



Collocation. First, we note that in view of Eq. (114), the present DSC-Galerkin statement reduces to a collocation 
one at the limit of a' 



lim 



N,Mi 



(x)0^i,,,(x)dx = <f (xO = 0, 



(119) 



where {x;} is the set of collocation points. However, in digital computations, we cannot take the above limits. It 
follows from the low pass filter property of the DSC trial functions, Eq. (115), that 



/ 

JQ. 



R^:^\x)0..':i{^)dx^R^:^\xi)^Q. 



(120) 



It can be proven that for an appropriate choice of S^, '^^ , the first approximation of Eq. (120) converges uniformly. 
The difference between the true DSC-coUocation, 



lim R^f{xi)^0, 



(121) 



and the Galerkin induced collocation, (120), diminishes to zero for appropriate DSC trial functions. 

Global and local. Global approximations to a function and its derivatives are realized typically by a set of truncated 
i^(a,6) function expansions. It is called global because the values of a function and its derivatives at a particular 
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point Xi in the coordinate space involve the full set of grid points in a computational domain fl. Whereas a local 
method does so by requiring only a few neighborhood points. In the present DSC approach, since the choices of M 
and/or M' are independent of N, one can choose M and/or M' so that a function and its derivatives at a particular 
point xi are approximated either by the full set of grid points in the computational domain or just by a few grid 
points in the neighborhood. In fact, this freedom for the selection of M endows the DSC algorithm with controllable 
accuracy for solving differential equations and the flexibility in handling complex geometries. 

Finite Difference. In the finite difference method, the differenti al op erator is approximated by difference operations. 
In the present approach, the DSC-coUocation expression of Eq. (12C) is equivalent to a 2M + 1 (or 2M) term finite 
difference method. This follows from the fact that the DSC approximation to the nth order derivative of a function 
can be rewritten as 



dx^ 



k+M 



(122) 



l=k-M 



where c^^ are a set of DSC weights for the finite difference approximation and are given by 



-kl.M 



d1 

dxi 



(123) 



Note that the standard finite difference scheme is obtained if 0* g..; is chosen as the Lagrange kernel as discussed 
earlier. Obviously, for each different choice of (t>^„, we have a different DSC-finite difference approximation. Hence, 
the present DSC approach is a generalized finite difference method. This DSC-finite difference was tested in previous 
studies |41|. When M — \, the DSC-finite difference approximation reaches its lowest order limit and the resulting 



matrix is tridiagonal. In this case, by appropriately choosing the parameter a, the present DSC weights c^^ can 
always be made exactly the same as those of the standard second order central difference scheme, i.e. -^^^ 0, for 
the first order derivative and -27,— for the second order derivative. Here A is the grid spacing. However, for 
a given numerical bandwidth, the DSC-finite difference approximation does not have to be the same as the standard 
finite difference scheme and can be optimized in a practical application by varying a. Another important choice of 
the DSC bandwidth is that M = iV, where N is the matrix length. Obviously, the computational matrix is no longer 
banded and this is a case we called a "global finite difference method" . These interesting features are illustrated by 
using a few numerical examples in the next section. 



IV. APPLICATION TO ELECTROMAGNETICS 



In this section, we examine the usefulness, test the accuracy and explore the limitation of the DSC approach for 
solving problems in electromagnetics. Many DSC kernels discussed in the previous sections are suitable for being used 
as DSC trial functions. For simplicity, we focus on three DSC kernels of the Dirichlet type, a regularized Shannon's 
kernel (RSK), 



M t ^ sinf(x-a;fe) 



a regularized Dirichlet kernel (RDK), 



sm 



Xk) 



(2L + l)sin 



IT X — Xk 



A 2L+1 



{x - Xfc)^ 

2a2 



exp 



{x - XkY 
2f72 



(124) 



(125) 



and a regularized Lagrange kernel (RLK) 



i=k+L 

n 

-k—L.i^k 



Xk 



■ exp 



{X - Xkf 



2(t2 



(126) 



for our numerical experiments. Note that the resolution is given by a = -Jj which is the frequency bound in the 
Fourier representation. The half bandwidth, M, can be chosen to interplay between the local limit and the global 
limit and is set to 40 in all calculations except for specified. Finally, L controls the order of the regularized Dirichlet 
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and Lagrange kernels and is set to 50 in calculations. It should be point out that, the selection oi L {L > M) is 
independent of the grid used in the present computation. 

Three different problems in computational electromagnetics, waveguide analysis, electrostatic analysis and 3D 
electromagnetic wave propagation, are selected to test the DSC algorithm. In all cases, we utilize the DSC algorithm 
to solve differential equations. Details of these computations are described in the next three subsections. Double 
precision is used in all calculations. 



A. Waveguide analysis 

EiGENMODE ANALYSIS The propagation of uniform plane waves is characterized by the transverse electromagnetic 
(TEM) nature of the wave. The electric and magnetic field intensities arc orthogonal to each other and to the 
direction of wave propagation. In a waveguide, the wave propagation is guided along given directions with particular 
characteristics determined by the structure of the wave guide. Although it is assumed that boundaries are perfectly 
conducting, the material within the waveguide is arbitrary and may include lossy or perfect dielectrics, conductor, 
etc. The properties of a waveguide with complex geometries can only be simulated numerically. To test the DSC 
algorithm for CEM applications, we first consider the problem of finding the eigenvalues that determine the parameters 
of waveguide modes, resonant frequencies of resonators, and many other physical parameters. Consider a rectangular 
waveguide, with propagation in the ^;-direction. For the transverse magnetic (TM) mode, the four field components 
Ex,Ey, Hx and Hy can be expressed in terms of Ez- In turn, Ez can be written as 

Ez{x,y,z,t) = E{x,y)e'^'''-'''\ (127) 

where E satisfies 



k^E = 0, (128) 



and vanishes on the boundary {E = 0) for the TM modes. The eigenvalue fc^ determines the phase parameter a 
through 

k"^ =u;^iye- a^, (129) 

where e and v are dielectric constant and magnetic permeability, respectively. 

To simplify the problem further, we take the computational domain as [0, IOtt] x [0, IOtt]. Three DSC kernels are 
used to compute the eigenvalue problem at three different grid sizes (36^, 24^ and 12^) so that the rate of convergence 
of the DSC approach for waveguide analysis can be examined. The computational bandwidth is set to the same as the 
grid size in this problem. When bandwidthes are M = 36, 24 and 12, the values of cr/A are chosen as 4.2, 3.2 and 2.65 
for the RSK and RDK and 2.95, 2.75 and 2.3 for the RLK. Results are obtained by a direct matrix diagonalization 
using a standard eigenvalue solver. Selected eigenvalues and absolute errors of the present numerical calculation are 
listed in TABLE I. 

We first note that the RSK and RDK have very similar performance in all computations. Hence, we only need to 
consider one of them in the rest of the computation. However, the performance of the RLK differs slightly from that 
of the RSK. 

For N = Nx = Ny = 12, there are only a maximum of 2.4 grid points per wavelength in each dimension. Obviously, 
if this is the case, it is impossible for any of the three kernels or any other method to catch up with the fast oscillations 
in high-order cigcnmodcs. All three kernels deliver excellent results when the mesh is refined to 24^ (about 4.8 grid 
points for each wavelength in each dimension). Machine precision is attained for the first 100 eigenmodes when the 
mesh is further refined to 36^. 

Field distribution In practical applications, the geometry of a waveguide can be very complex so as to achieve 
certain required properties. Therefore, it is important for a computational algorithm to be fiexible enough in handling 
complex boundary conditions. In this study, we consider two waveguides with geometries having a T-shape and an 
E-shape. These boundaries are designed to test the capability of the DSC algorithm for complex waveguide simulation. 
We are interested the field distribution of TM (and/or TE) modes of these waveguides. The RSK is used for both 
cases and the value of cr/A is set to 4.8. A total of 50 DSC-coUocation points are used in each dimension. Field 
distributions across the cross section for eigenmodes 1, 2, 3 and 4 are shown in Fig. la. These eigenmodes are quite 
localized and do not give the full shape of waveguide confining geometry. However, as plotted in Fig. lb, higher-order 
eigenmodes, 17, 18, 19 and 20, clearly reflect the waveguide geometry. Similar features are observed in the E-shape 
waveguide as shown in Fig. 2. The first four modes of E-shape waveguide are given in Fig. 2a. Modes 17, 18, 19 and 
20 are depicted in Fig. 2b. 
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B. Electrostatics analysis 



Potential function provides a useful alternative representation of electromagnetic fields. In many situations, it is 
more convenient to use an auxiliary function to describe field properties rather than to use full electric and magnetic 
field variables. In particular, a scalar potential can be used in a source- free domain and/or in domains which contain 
scalar sources such as charge densities. 

The electric scalar potential is one of the most widely used potential functions. If the rate of change of the magnetic 
flux vanishes, the electric field is irrotational 

V X E = 0. (130) 
In such a case, the electric field can be written as the gradient of a scalar potential 

E = -VV. (131) 



For static fields, the divergence of Eq. (131) gives 



VV = --, (132) 

e 



where p is the volume electric charge density and e is the electric permittivity. Here, the Maxwell equation V • E 



_ £ 



is used for the derivation. Equation (132) can be used as a starting point for many useful theoretical descriptions. 

To demonstrate the DSC algorithm for electrostatic calculations, we first consider a simple problem in two spatial 
dimensions. The geometry of the problem consists of a conducting square box of Im x Im in the x — y plane. The 
box is infinitely long in the z — direction and thus, this can be regarded as a two-dimensional problem. The boundary 
condition of the problem is that the top side is isolated from all other sides with a potential of V{x, 1) = lOV^. All 
other three sides are connected to the zero potential (y(0,j/) = V{l,y) ~ ^(2:, 0) = 0). The material inside the box 
is the free space with permittivity £0 = x 10"^. 

We next consider an extension to the first problem. The same box with the same boundary condition is used except 
that there is a 0.18m x 0.16m area of charge density inside the box. The charge density is uniformly distributed 
within a material of relative permittivity Sr — 100 and has a value of 1.0 x 10~^[C/m^]. 

In both cases, we choose 32 points {N ~ 32) in each dimension (Aa; = Ay — 0.0322to). The DSC-collocation 
method is very convenient for handling these potential function problems. The solution is obtained by solving the 
coupled algebraic equations (implicit scheme) with all boundary conditions being implemented in the operator part 
(£) of the coupled equations. The charge density is treated as a source term (/) in the coupled equations. For the 
first case, there is no charge density in the inner region and the Laplace equation is directly solved. The results are 
plotted in Fig. 3. 



For the second case, Eq. (132) is modified as 



VV(x,t/) = -^. (133) 

The charge density at the area of x : 0.41to — 0.59m, y : 0.72m — 0.88m and the boundary conditions are directly 
imposed in a set of coupled algebraic equations, which are generated by using the DSC-collocation formulation. A 
standard linear algebraic equation solver is used to attain results as plotted in Fig. 4. 

C. Electromagnetic wave propagation 

In an uncharged homogeneous media, propagation of electromagnetic waves is governed by the wave equation 

d^W 1 



-^'W = 0, (134) 

where PF(r,t) can be either the electric field or the magnetic field. To illustrate the potential of the DSC algorithm 
for electromagnetic wave simulations, we consider the following initial value problem 

W{x, y, z, 0) = sm{axX + ayy + a^z), (135) 

with periodic boundary conditions in all directions. This problem has an analytical solution given by 
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W{x, y, z, t) — sin 



' -(- Q;2 _|_ Q,2 



(136) 



To further simplify notations, we set Ux = cty — oiz — — ^- The initial wave is propagated in a cubic size (IOtt)'^. 
The fourth order explicit Runge-Kutta scheme is used for time discretization. Sufficiently small time increments (Ai) 
are used so that the major errors are d ue t o the spatial discretization. The RSK and RL K ar e used in the finite 



difference manner, i.e. according to Eq. (122), and DSC weights are calculated by using Eq. (|123| ) once for the whole 
computation. We compute the problem by using a variety of grid points N [N = Nx = Ny = Nz) and computational 
bandwidth M (M = Mx = My = Mz) values to test the accuracy and rate of convergence of the DSC algorithm and 
to demonstrate the unified features of global and local methods in the DSC algorithm. When M =36, 24, 12, and 6 
the values of a/ A are 4.2, 3.2, 2.65 and 2.0 for the RSK and 2.95 2.75, 2.3 and 1.85 for the RLK. The Loo errors at a 
number of time units are listed in TABLE II. 

When N = 24, the number of grid points per wavelength is 4.8 in each dimension (Note that typically the Yee 
algorithm uses about 18 points per wavelength). This is a typical case of under sampling and it is very difficult to 
achieve high computational accuracy by any means. We choose three different M values. For M — 6, results from 
the two DSC kernels are accurate enough. The accuracy improves as the M value is increased from 6 to 12 and 24. 
The last case, M = N = 24, corresponds to a global treatment (i.e., a "global finite difference") and its results are 
significantly better. 

In the case oi N = 36, it is on an average about 7.2 grid points per wavelength in each dimension. The computational 
accuracy increases as the computational bandwidth increases. It is noted that at M = 12 and M = 24, the accuracy 
does not improve much from the previous computations {N — 24; M = 12, 24) and is obviously limited by the 
computational bandwidth. However, it is interesting to note that the DSC algorithm attains extremely high accuracy 
at the global limit of Af = iV = 36. To the best of our knowledge, these are the highest precision results obtained for 
this problem, so far. 



V. CONCLUSION 



In conclusion, this paper introduces the discrete singular convolution (DSC) algorithm for computational electro- 
magnetics (CEM). The computational philosophy of the DSC algorithm is studied. Many sequences of approximations 
to the delta distribution, the "universal reproducing kernel" , are constructed either as bandlimited reproducing ker- 
nels or as approximate reproducing kernels. A rcgularization procedure based on the distribution theory is utilized 
to improve the regularity and localization of the DSC kernels. Systematic treatments of derivatives and boundary 
conditions are proposed, as they are required for a computational algorithm. 

We explore the unified features of DSC algorithm for solving differential equations in the framework of the method 
of weighted residuals. It is found that several computational methods, such as global, local, Galerkin, collocation and 
finite difference methods, can be derived from a single starting point by using DSC trial functions and test functions. 
The unification of local and global methods is realized via the DSC approach. It is well known that accuracy is 
crucial to many scientific computations, such as simulation of turbulence and high-frequency analysis of radar cross- 
section. Whereas, small-matrix-band approximations are vital to large scale computations. The DSC approach 
provides a controllable numerical accuracy by an appropriate selection of the matrix bandwidth. Therefore, by using 
the DSC algorithm, the computational efficiency can be easily optimized against accuracy, convergence, stability and 
simplicity. A collocation algorithm is deduced from a Galerkin statement and is called a Galerkin-induced collocation. 
The equivalence of Galerkin and collocation enables us to evaluate these two conventional numerical algorithms on an 
equal footing. The connection is made between finite difference and other methods, such as collocation and Galerkin. 
The DSC algorithm can be regarded as a generalized finite difference method, which becomes a "global finite difference 
method" by an appropriate choice of the computational bandwidth. These results can be used as a guide for the 
selection of computational methods and for the design of numerical solvers for practical applications. 

The potential of the DSC algorithm for computational electromagnetics is explored by performing a number of 
numerical experiments. The first example is about waveguide eigenmode analysis, an eigenvalue problem in a square 
domain of (IOtt)^. This problem has an analytical solution and hence can be used for testing the accuracy of potential 
numerical methods. We formulate the problem in the DSC-coUocation approach and attain results by using three 
sets of grid points, 12^, 24^ and 36^, to examine the speed of convergence of the DSC algorithm. Three typical DSC 
kernels, a regularized Shannon's kernel, a regularized Dirichlet kernel, and a regularized Lagrange kernel are utilized 
to illustrate the usefulness of the proposed algorithm. It is found that the accuracy improves dramatically as the 
number of grid point is increased. In particular, it reaches the machine precision for the first 100 cigcnmodes when 
the mesh is refined to 36^. Having built enough confidence for eigenmode computations, we have tested the ability of 



24 



the DSC algorithm in handhng complex boundary conditions. In this regard, we have used the DSC algorithm for TM 
field distributions of T-shape and E-shape waveguides. We found that the first few eigenmodes are very localized due 
to complexity in geometry. However, higher-order eigenmodes exhibit global characteristics. The confining geometry 
of the waveguide can be seen from mode patterns. 

In the second example, we consider potential function analysis of electrostatics in two spatial dimensions. Com- 
putationally, this is a non-uniform boundary value problem. Only one side of the four boundaries has a non-zero 
potential value. We have treated the problem via the DSC-coUocation approach with a mesh of 32^. Since the DSC 
algorithm is a local method in general, boundary conditions and charge density are easily implemented in a set of 
coupled algebraic equations. The DSC approach has proven to be very efficient and robust for this problem. 

In the last example, we have tested the DSC algorithm for three-dimensional electromagnetic wave propagations. 
The wave equation is integrated over a long time with periodic boundary conditions. Since the problem admits an 
analytical traveling wave solution, the performance of the DSC algorithm can be examined objectively. This problem 
is rather sensitive to the computational accuracy and stability. The temporal discretization is carried out by using the 
standard fourth order explicit Runge-Kutta scheme. The DSC algorithm is utilized for spatial discretization in the 
domain of (IOtt)^. The unified feature of the DSC algorithm for both global and local approximations is illustrated 
in this problem by varying the computational bandwidth AI for a given matrix size N. In this calculation, the DSC 
algorithm is used as a generalized finite difference method. In particular, when the bandwidth equals the matrix length 
(M = N), the DSC algorithm gives rise to an interesting "global finite difference" treatment. The DSC algorithm 
can provide controllable accuracy by varying the parameter Af , hence, one of the advantages of the DSC algorithm 
is its robustness, and allows the selection of desirable accuracy for a given problem without any need to change one's 
computer code. All numerical solutions are found to be stable over a long time integration. A sharp improvement 
in the numerical accuracy is observed when the mesh is refined from 24^ to 36^, i.e. from 4.8 points per wavelength 
to 7.2 points per wavelength in each dimension. It is found that the accuracy of the DSC algorithm is at least 12 
significant figures up to 22 time units when the computational bandwidth is relatively high (Af = N = 36). This 
result is consistent with our previous theoretical error estimation |Q. Present work indicates that the DSC algorithm 
is a promising and potential approach for computational electromagnetics. 
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Figure Caption 

FIG. la. Contour and mesh plots of the first four eigenmodcs of the T-shape waveguide. 
FIG. lb. Contour and mesh plots of eigenmodes 17, 18, 19 and 20 of the T-shape waveguide. 
FIG. 2a. Contour and mesh plots of the first four eigenmodes of the E-shape waveguide. 
FIG. 2b. Contour and mesh plots of eigenmodes 17, 18, 19 and 20 of the E-shape waveguide. 
FIG. 3. Contour and mesh plots of the potential field distribution in a square box with non-uniform boundaries. 
FIG. 4. Contour and mesh plots of potential field distribution in a square box with non-uniform boundaries and 
an area of charge density. 
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TABLE I. Errors for waveguide modes 



N 


Eigenmode 


RSK 


RDK 


RLK 


12 


10 


8.4(-8) 


8.8(-8) 


2.1(-7) 




20 


5.2(-7) 


5.1(-7) 


4.7(-7) 




30 


1.9(-2) 


1.9(-2) 


2.0(-2) 




40 


l.l(-l) 


l.l(-l) 


l.l(-l) 




50 


1.5(-1) 


1.5(-1) 


1.5(-1) 




60 


1.9(-1) 


1.9(-1) 


1.9(-1) 




70 


2.7(-l) 


2.7(-l) 


2.7(-l) 




80 


3.2(-l) 


3.2(-l) 


3.2(-l) 




90 


3.7(-l) 


3.7(-l) 


3.7(-l) 




100 


4.2(-l) 


4.2(-l) 


4.2(-l) 


24 


10 


1.8(-14) 


3.1(-15) 


9.7(-13) 




20 


9.5(-15) 


1.6(-14) 


3.2(-12) 




30 


8.5(-13) 


8.2(-13) 


1.4(-12) 




40 


8.7(-13) 


8.5(-13) 


5.1(-13) 




50 


3.8(-10) 


3.7(-10) 


1.3(-11) 




60 


1.4(-8) 


1.3(-8) 


3.5(-10) 




70 


8.0(-8) 


7.8(-8) 


6.3(-9) 




80 


8.0(-8) 


7.8(-8) 


6.3(-9) 




90 


7.6(-10) 


7.5(-10) 


3.0(-ll) 




100 


6.4(-9) 


6.3(-9) 


3.6(-10) 


36 


10 


1.2(-14) 


2.1(-14) 


2.9(-14) 




20 


4.3(-14) 


2.1(-14) 


1.2(-15) 




30 


1.4(-14) 


2.8(-14) 


7.7(-15) 




40 


2.2(-16) 


3.9(-14) 


5.0(-15) 




OU 










60 


7.4(-14) 


2.0(-14) 


3.8(-14) 




70 


2.4(-14) 


7.3(-15) 


2.9(-14) 




80 


2.2(-14) 


4.2(-14) 


6.7(-15) 




90 


1.3(-14) 


6.4(-15) 


4.4(-14) 




100 


3.7(-14) 


1.8(-15) 


3.3(-14) 
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TABLE II. Loo errors for solving the 3D wave equation 



N 


M 


At 


t 


RSK 


RLK 


24 


6 


0.5 


1 


6.2(-4) 


5.1(-4) 








4 


l.l(-3) 


1.5(-3) 








7 


1.7(-3) 


2.4(-3) 








10 


2.3(-3) 


3.3(-3) 








13 


2.8(-3) 


4.3(-3) 








16 


3.5(-3) 


5.2(-3) 








19 


4.1(-3) 


6.1(-3) 








22 


4.8(-3) 


7.0(-3) 




12 


0.2 


1 


1.2(-5) 


1.2(-5) 








4 


2.8(-5) 


3.1(-5) 








7 


4.5(-5) 


5.1(-5) 








10 


6.1(-5) 


7.1(-5) 








13 


7.8(-5) 


9.0(-5) 








16 


9.4(-5) 


1.1 (-4) 








19 


l.l(-4) 


1.3(-4) 








22 


1.3(-4) 


1.5(-4) 




24 


0.05 


1 


4.3(-8) 


5.0(-8) 








4 


7.6(-8) 


2.0(-7) 








7 


l.l(-7) 


3.4(-7) 








10 


1.4(-7) 


4.9(-7) 








13 


1.8(-7) 


6.4(-7) 








16 


2.1(-7) 


7.8(-7) 








19 


2.4(-7) 


9.3(-7) 








22 


2.7(-7) 


l.l(-6) 


36 


12 


0.2 


1 


1.3(-5) 


1.3(-5) 








4 


4.8(-5) 


4.7(-5) 








7 


8.6(-5) 


8.4(-5) 








10 


1.2(-4) 


1.2(-4) 








13 


1.6(-4) 


1.5(-4) 








16 


2.0(-4) 


1.9(-4) 








19 


2.3(-4) 


2.2(-4) 








22 


2.6(-4) 


2.5(-4) 




24 


0.05 


1 


5.2(-8) 


5.2(-8) 








4 


2.1(-7) 


2.1(-7) 








7 


3.6(-7) 


3.6(-7) 








10 


5.2(-7) 


5.2(-7) 








13 


6.8(-7) 


6.8(-7) 








16 


8.2(-7) 


8.2(-7) 








19 


9.8(-7) 


9.8(-7) 








22 


l.l(-6) 


1.1 (-6) 




36 


0.001 


1 


1.3(-14) 


1.4(-14) 








4 


3.0(-13) 


3.0(-13) 








7 


i .0\-Lo ) 


( .ol^-io ) 








10 


3.5(-14) 


2.7(-14) 








13 


1.7(-12) 


1.7(-12) 








16 


3.3(-12) 


3.3(-12) 








19 


3.9(-13) 


4.1(-13) 








22 


4.1(-12) 


4.1(-12) 
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